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Abstract 

In this paper we extend previous work deriving dynamic equations governing infectious disease spread on networks. The 
previous work has implicitly assumed that the disease is initialized by an infinitesimally small proportion of the population. 
Our modifications allow us to account for an arbitrarily large initial proportion infected. This helps resolve an apparent 
paradox in earlier work whereby the number of susceptible individuals could increase if too many individuals were initially 
infected. It also helps explain an apparent small deviation that has been observed between simulation and theory. An 
advantage of this modification is that it allows us to account for changes in the structure or behavior of the population 
during the epidemic. 
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Introduction 

The mathematical study of infectious disease spread has 
contributed significantly to our ability to design effective interven- 
tions to reduce disease spread. Most of the earliest models were 
based on the assumption that disease transmission occurs as a 
Poisson process and each transmission reaches an individual 
chosen randomly from the population. This implicitly assumes that 
partnership duration is very brief. These models have been 
modified to account for a number of different effects, such as 
demographic groups [1]. 

More recendy, attempts have been made to incorporate the 
"network" structure of the population (see, e.g., [2]). Typically 
these focus on trying to understand the role played by "high- 
degree" individuals (those individuals with many contacts), and 
they come in one of two flavors: they either continue the 
assumption of fleeting partnerships (the timescale of individual 
transmissions is long compared to the timescale of individual 
partnerships) [1,3-6], or they take the opposite limit in which the 
partnership network is static (the timescale of the epidemic is short 
compared to the timescale of individual partnerships) [7-14]. 
These two approaches do not address the intermediate regimes. 

Recent work has shown that for susceptible-infectious-recovered 
(SIR) models, it is possible to unify these two approaches with an 
"edge-based compartmental model" (EBCM) that allows partner- 
ship duration to range continuously from zero to infinite [15-17] 
[for susceptible-infectious-susceptible (SIS) models, the picture is 
more complicated, see for example [18]]. The resulting models are 



low-dimensional and contain many standard models as special 
cases [17]. Unfortunately, these models are derived under the 
assumption that the initial proportion infected is infinitesimally 
small (while the absolute number infected is sufficiendy large that 
the dynamics are deterministic). It is assumed that by the time the 
equations are used, any early transients have died away. A 
consequence of this assumption is that the models break down if 
TZq < 1 (that is, if the average number of infections caused by an 
infected individual early in the epidemic is less than 1) or if the 
initial proportion infected is not negligible. 

The failure if the initial proportion infected is not negligible was 
observed by [14]. This paper used an early (static network) version 
of the equations of [15] from [7] and compared them with 
simulation. A small discrepancy in final sizes was noted. This 
discrepancy was not present for equations of [19], a system 
requiring O(M) equations where M is the maximum degree, or 
for another system introduced in [14] which required 0(M 2 ) 
equations. 

In the remainder of this paper, we generalize the EBCM 
equations for the spread of infectious disease assuming a finite 
proportion of the population is initially infected. We test the 
resulting equations against simulations, analyze their predictions 
for different disease scenarios, and investigate the cause of the 
discrepancies found in previous work. For simplicity we focus on 
the static network limit. The method we introduce is straightfor- 
ward to adapt to dynamic networks. 
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Analysis 

We modify the approach of [15] which assumed an infinitesimal 
initial proportion infected. We adapt the approach to consider a 
wide range of possible initial conditions. We assume that the 
dynamics of the epidemic may be treated as deterministic, which 
means we assume the population is very large and the initial 
number infected is large enough for the epidemic to behave 
deterministically. The assumption that behavior is deterministic 
may be understood qualitatively as equivalent to the claim that no 
individual has a large enough effect to alter the dynamics of the 
disease at the population-scale: the time (or even if) a single given 
individual becomes infected has a negligible impact on the 
proportion infected. If stochastic effects are still important but 
TZo > 1, then these equations may become accurate at a later time 
once sufficient numbers are infected. 

We assume the population consists of JV» 1 individuals. Each is 
assigned a degree k (independently of degrees of other individuals) 
with probability P(k) where P defines a probability distribution on 
the non-negative integers. The network is wired together using the 
"Configuration Model" (or "Molloy-Reed") approach [11,20]: 
each individual is assigned a number of stubs (or half-edges) equal 
to its degree. Pairs of stubs are then wired together to form edges/ 
partnerships. It is likely that this algorithm produces a handful of 
self-loops or repeated edges, but although they may be present, 
their density [i.e., the probability a given individual is involved in a 
self-loop or repeated edge) goes to zero like l/N. 

We define a test individual u to be a random individual chosen at 
time 0. Because we assume that the spread is deterministic, this 
means that the probability a is in a given state is equal to the 
proportion of the population in that state. So we focus on 
calculating the probability u is susceptible, infected, or recovered. 
We modify u so that it does not transmit to any of its partners if 
ever infected. This assumption does not affect the probability u is 
in any given state, but it does prevent a correlation between the 
statuses of different partners which would be caused by infection 
traveling through u. This allows us to treat the partners of u as 
independent and so each partner of u may independently transmit 
to u. It is important to note that this assumption has no impact on 
the probability u is in any given state and therefore, it does not 
affect our calculation of the proportion of the population in each 
state. Further discussion of the test individual is in [21]. 

Variables and Parameters 

We introduce our variables, our parameters, and their 
definitions in table 2. The starting point is the test individual u. 
The remaining variables and parameters can be broadly divided 
into four groups. 

• S, I, and R denote the proportion of the population in each 
state, or equivalendy the probability that the test individual u is 
in each state. 

• 9, (j) s , and (j) R give the probability a partner of u has a given 
status and the probability the partner has transmitted to u: 9 is 
the probability the partner has not transmitted to u and 0s, 0/, 
and (j) R give the probability the partner has not transmitted to 
u and is susceptible, infected, or recovered respectively. 

• P(k) tells us the probability a random individual has degree k, 
while S(k,0) tells us the probability a random individual of 
degree k is initially susceptible. The function \jj(x) = S(k,Q) 
P(k)x k encodes P and S. We define (K} = Y] k kP(k) t0 be 
the average degree. 

• We have two disease parameters to consider: /?, the 
transmission rate, and y, the recovery rate. 



Given our definitions, ij/(9(t)) is the probability u is susceptible 
at time t. By noting that \jj(9(t)) = S(t), we will be able to close our 
system of equations. 

The main distinction between this approach and the previous 
EBCM approach [15] is that we use just the initially susceptible 
individuals to define \jj while the earlier work took 
{ ! J (. x ) = *l2iP(k) xl< , the probability generating function for the 
degree distribution. The earlier work then assumed the disease had 
already been spreading prior to time 0 and defined 9 to be the 
probability that a partner has never transmitted [so 8(0) < 1] 
whereas here we take 9 to be the probability that a partner has not 
transmitted given that it had not prior to time 0 [so 9(0)= 1]. 

Equation Derivation 

We will find a closed system of equations based on these 
variables. We begin by looking at S(t). If the test individual u has 
degree k and is susceptible at t = 0, then the probability it is 
susceptible at some later time is 9(t) k . If we do not know k or 
whether u is susceptible at ? = 0, then the probability u is 
susceptible at time t is the sum over all k of the product of the 
probability u is initially susceptible S(kft) with the probability u is 
still susceptible 9 k . We have S(t)= J2 k S(k,0)P(k)9(t) k = \j/(9(t)). 
Thus we conclude 

S(t) = ^(9(t)) 

We know that R solves R = yl. We also have a conservation rule 
that S + 1 + R=\, so 7=1— S — R. Thus our equations are 



I=l-S-R 



R = yl 

Assuming 9(t) is known, then this system with an initial 
condition for R completely defines S, I, and R. This is shown in 
the flow diagram in figure 1 . Other formulations are possible, for 
example S = \J/(9), 1= — \j/'(9)9 — yl, R = yl. However we find 
our system to be preferable because it minimizes the number of 
differential equations. 

In order to close this system of equations we need an equation 
giving 9. Recall that 9(t) is the probability a partner v of u that had 
not yet transmitted to u by time 0 has still not transmitted by time 
/. This is broken into three disjoint sub-compartments 



S = 




I 




R 







Figure 1. Flow diagram showing the flux of individuals 
between the different compartments. Because we have an explicit 
expression for S, if we know 0 we do not need to explicitly determine 
the flux from S to /. 
doi:1 0.1 371/joumal.pone.01 01 421 .g001 
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8 = <j>s J r$i J r<j> R based on whether v has not transmitted and is 
susceptible, infected, or recovered. Because ^/ is the probability v 
has not transmitted and is infected, it is straightforward to see that 
8 = — Pfij. So if we can find in terms of 8, then we arrive at a 
single equation for 8, which can be used to provide 8 for the S, I, 
and R equations. 

To do this, we use the fact that <l>i = 8 — (j) s — (l> R and find <j> s and 
(fi R in terms of 8. We turn to figure 2. The recovery rate is y and 
the transmission rate is /J, so we have (j) R = — yd/ ft. We can 
integrate this, and using the fact that 8(0) =1 we find 
<j) R = y(\—8)/fl + (l> R (0). To find tj> s in terms of 8, we note that 
the probability u has an edge to a node that is susceptible at time 
t = 0 is 0j(O). The probability the susceptible partner has degree k 
is kP{k)S(k,0)/ J2kP(k)S(kfi), so the probability an initially 
susceptible partner is susceptible at some later time is 
J2 k k P(k)S(k,0)8 k - 1 /J2k k P(k)S(k,0) = \j/'(8)/\jj'(\). Thus 
0 S (^ = 0 S (O)V'(0(O)/>A'(1)- w e arrive at 



^ = 0-*s(O)^-|(l-fl)-^(O) 



and 8= — P(j)j becomes 



>+/ty S (O)^+Y(l-0)+/ty*(O) 



with 8(0) = 1 . This completes our system. 
Our final closed system of equations is 



>+/ty S (O)^+Y(l-0)+^ R (O) 



R = yl, S = ljf(9), I=l-S-R 



(2) 



where 8(0) = 1, and R(0) is given by the initial conditions. These 
equations lead to earlier equations of [7,15,22] if TZo>\ and 
1 —<j> s (0), </>j(0), <I> R (0) and 1 — 0(0) are all infinitesimally small. If 
72-0 < 1, the error caused by the approximations <^(0) = 0 and 
(fi s (0)= 1 is comparable to the actual number infected. 

Final size relation. The final size relation assuming small 
initial condition is well-known [11]. The final size relation for 




1 - 8 



Figure 2. Flow diagram for the flux of partners through 
different states. The top three boxes § s , <j>j, and <j> R make up 0 and 
represent the different states the partner can be in if it has not 
transmitted. The lower box 1—0 is the probability the partner has 
transmitted. 

doi:10.1371/journal.pone.0101421.g002 



larger initial conditions has recently been found [21] in a more 
general case not assuming constant transmission and recovery 
rates. It can be derived easily for this model by setting 8 = 0. We 
find 



0(oo)= 



s (0) 



not*)) 
go) 



MO) 



+y 



(3) 



J R(oo)=l-l>(0(oo)) 



(4) 



Model Validation 

We now compare our model with simulations for populations 
that satisfy the Configuration Model/Molloy-Reed model as- 
sumptions. Although an earlier version of these equations was 
found to have minor discrepancies [14], we show that once we 
appropriately account for the initial condition, the calculation 
becomes correct. 

Final Size Comparison. To show that our new equations 
accurately calculate the impact of the initial conditions, we first 
consider epidemic spread in networks with the same degree 
distribution as in [14] (table 1), but with varying numbers initially 
infected and varying population sizes. We then consider the 
impact of selecting high or low degree nodes as the earliest infected 
individuals, using networks whose degree distributions more 
clearly show the impact of biased selection of the initial 
individuals. 

We run a large number of simulations for each number of initial 
infections. For each simulation we generate a new network. Our 
simulation technique is similar to those recendy described by 
[8,19,23]. In the Configuration Model framework, each node is 
assigned a degree, nodes are given stubs (or half-edges), and then 
stubs are randomly paired together. In the simulations we use, 
each node is assigned a degree, nodes are given stubs, and then the 
disease begins to spread in the network before stubs are paired. 
Each time the disease transmits along a stub that stub is joined to a 
randomly selected as-yet-unpaired stub. If the partner is suscep- 
tible, then it becomes infected. If not, nothing happens. Once stubs 
are paired they remain in their edge. This approach is equivalent 
to constructing the network in advance and then following the 
disease, but it is more efficient computationally because it only 
constructs those parts of the network the disease traces. 

Randomly selected initial infections. We first consider 
varying numbers of randomly chosen initial infected individuals. 
In figure 3 we take the degree distribution from table 1. 

We randomly select a proportion p of the population to initially 
infect. We have S(k,0) = \—p for all k, so >p(x) = ( 1 — p) 
^ k P(k)x k . Similarly we have (/) s (0)=\—p. Because the 
epidemic begins with no recovered individuals, we take <I> R (0) 
IsO. We take ft = 0.1 and y = 0.2 (though all that matters for the 
final size is their ratio). 

We take populations of 100, 1000, and 10000 and perform 
many simulations. We compare the final sizes observed with the 
final size relation of equations (3) and (4). The equations are 
derived in the infinite population limit, but in figure 3 we see that 
even with populations of only 100 they give a good prediction of 
the observed behavior. As the population size increases, the 
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Table 1. The degree distribution used in simulations in [14]. 







k 


P(k) 


1 18.118 x 10- 3 


2 


72.536 x 10-- 1 


3 


145.222 x 10-- 1 


4 


194.589 x 10- 3 


5 


195.962 x 10- 3 


6 


156.857 x 10- 3 


7 


105.280 x 10- 3 


8 59.713 xl0~ 3 


9 


30.066 x 10- 3 


10 


21.657 x 10" 3 





doi:l 0.1 371 /journal.pone.Ol 01 421 .t001 



distribution becomes narrower and the simulations collapse more 
tightly around the prediction. 

Biased initial infections. To show that the approach we 
have derived can also be applied to cases where the initial infected 
individuals are selectively chosen based on their degree, we use a 
different degree distribution that helps highlight the effect. We take 
P(\) = P(9)= 1/2. We consider two options. In the first approach, 
individuals with higher degree are preferentially selected. To do 
the selection, we choose an individual with probability propor- 
tional to the square if its degree, and infect it. We repeat this until 
a proportion p of the population is infected. In the second 
approach individuals are chosen with probability proportional to 
the square of their inverse degree until a proportion p is infected. 
We take 0 = 0.1 and y = 0.6. 



Using these rules, we clearly see that S(k,0) is not uniform. Instead, 
for the case where individuals are selected with probability 
proportional to their squared degree, we find that S(k,Q) = a* 2 where 

a solves J2k p ( k )^ = 1 ~P- We fmd ta<P) = E kS(k,0)P(k)o^ / 
^P(k). In the case where individuals are selected with probability 
inversely proportional to their squared degree, we find that 
S(kfi) = a i l kl where a solves ^PQ^a 1 ^ = l-p, and <j> s (0) = 
£ k S(k,0)P(ky/ k2 1 Y, k k P(k). 

We compare predictions and simulations in populations of 1000 
individuals in figure 4. In the limit of a negligible initial proportion 
infected, the final size of epidemics in these networks is about 4%. 
As we increase the number of initially infected individuals, we 
increase the final size because of these individuals and because of 
the additional infections they lead to. At small amounts, increasing 
the initial number of high degree nodes has a much larger impact 
on the final size because they cause more additional infections. 
However, as the amount of infection initially present is increased 
this effect becomes less important: the high degree individuals 
would become infected anyway. So the largest gain in final size 
comes from infecting low degree individuals who would not 
receive an infection from their partners. The "kinks" that occur 
just above 50% initially infected are because effectively all 
individuals of high (left) or low (right) degree are initially infected. 

Dynamic Calculation. We now look at the performance of 
the dynamic equations. The dynamic prediction is more easily 
affected by noise than the final size prediction, so we use larger 
population sizes. We again take the degree distribution of table 1. 
We begin with 5% infected, either randomly chosen, or chosen as 
before proportional to the square of the degree. A comparison of 
individual simulations with theoretical predictions is in figure 5. 
The theory accurately predicts the dynamics of epidemics in large 
populations, but there is significant stochasticity in smaller 



Table 2. The variables we need to calculate the epidemic 
population and modified so that it cannot infect others, 


dynamics. In all of these u is a test individual: randomly chosen from the 
although it can become infected. 




Variable/parameter 


Definition 


Test Individual u 


A randomly member of the population chosen at time t = 0 who is prevented 
from transmitting to its partners. 


Sit) 


The proportion of the entire population that is susceptible. 


m 


The proportion of the entire population that is infected. 


R{t) 


The proportion of the entire population that is recovered. 


0(t) 


The probability a random partner v of u that did not transmit to u by time 0 has 
not transmitted to u by time /. 


hit) 


The probability a random partner v that did not transmit to u by time 0 is 
susceptible at time t. 


M0 


The probability a random partner v that did not transmit to u by time 0 is 
infected at time t but has not transmitted to u. 


Mr) 


The probability a random partner v that did not transmit to u by time 0 is 
recovered at time t and never transmitted to u. 


m 


The probability an individual has degree k. 


<K)=Y, k kP(k) 


The average degree. 


S(k,0) 


The probability an individual with degree k is initially susceptible. 


H9(t))=j: k S(k,0)P(k)9(t) k 


The probability that the test individual u is susceptible at time t. In a large 
population this should equal S{t). 


ft the per-edge transmission rate 


y 


the per-individual recovery rate 


doi:1 0.1 371 /journal.pone.Ol 01 421 .t002 
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0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 

f> P P 

Figure 3. Results of simulations for N = 100, 1000, and 10000 individuals. The solid curve gives our prediction for the final sizes of epidemic in 
a large population. Colors are log scale giving probability of that particular epidemic size. Each simulation is for a new network generated using the 
P(k) from table 1, with /? = 0 . 1 and y = 0.2. We randomly select a proportion p of the population to initially infect and compare final size with the 
prediction of theory. The number of simulations for each p for N= 100, 1000, and 10000 was 50000, 11500, and 3000 respectively. To show that this is 
sufficient to resolve the distribution, for .175<p<.225 there were 2000000, 50000 and 10000 simulations performed respectively for each N. This 
only slightly improves the tails of the distribution. Note that when the initial number (not proportion) of infections is small, a large fraction of 
simulations end without an epidemic. 
doi:10.1371/journal.pone.0l0l421.g003 



populations. The theory is able to capture the impact of the biased 
selection of initial infections. 

Intervention Impact 

We can use our equations to compare the impact of several 
interventions. We consider an epidemic spreading in the 
population, and at some intermediate time we introduce a change 
in the disease or population. Because the system changes at a time 
with a non-negligible amount of infection in the population, the 
equations derived assuming a negligible proportion infected fail. 

Consider a population in which P(4) = P(5) = P(6) = 1/3. 
Assume we initially infect a small, randomly chosen proportion 
of the population, p at t = 0. Thus we have 
iA(.y) = (1-p)(a- 4 + x 5 +x 6 )/3, ( j> s (0)=\-p, <t> I (G) = p, and 
<^(0) = 0. We take /?=1 and y= 1/2. 

We consider three interventions that may be introduced at time 
/ 1 . All are aimed at "halving" the transmission rate, but they do 
this in different ways. In mass-action based models, these would all 
have the same effect. We can clearly identify differences using our 
approach. 

• An intervention that reduces ft by a factor of 2. 

• An intervention that reduces /? so that per-contact transmission 
probability P/(P + y) is reduced by a factor of 2. 




P 



• An intervention that eliminates half of the partnerships 
randomly. 

The distinction between the first two comes from the fact that 
partnerships have duration, so reducing (I by a factor of 2 does not 
reduce the infection probability by a factor of 2. The expected 
number of transmissions an individual sends to a given partner is 
Ply, but once the partner is infected, the subsequent transmissions 
are irrelevant. If we use mass action assumptions however, each 
transmission is to a replacement partner. In each case, halving /? 
halves the total number of transmissions, but when partnerships 
have nonzero duration, a larger proportion of transmissions have 
no effect. The probability of transmitting at least once along a 
static partnership is /?/(/? + y). So to reduce infection probability 
by a given factor requires a larger reduction to /}. Note that the 
work of [24] suggests that in Configuration Model networks the 
final size of our second and third intervention will be the same, 
(but that in clustered networks it will be different). 

We will demonstrate our approach in all three cases, restarting 
the calculations when the intervention is put into place. In all 
cases, this allows us to use the conditions at 1 1 to predict the final 
size. We take \l/ 0 (x) = \j/(x), 6q, (j> s $, 0/ jO ; an d <^j?,o to correspond 
to time less than t \ . We use a subscript of 1 for times after t \ . We 




Figure 4. Epidemic final sizes in population of 1000 individuals with half having degree 9 and half with degree 1. The disease 
parameters are /3 = 0.1, y = 0.6. Results of simulations having initial infections chosen with probability proportional to square of degree (left) or 
inverse square of degree (right). For each initial number of infections, 22500 simulations were performed, each with a different network. For the range 
0.175 <p<0. 225, 200000 simulations were performed to give insight into how well resolved the distribution is. Note that for small numbers of initial 
infections, epidemics are less likely when the lower degree individuals are chosen. 
doi:10.1371/journal.pone.0101421.g004 
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0.20 



0.15 



0.05 



0.00 





Figure 5. A comparison of the observed and predicted number of infections from simulations. Left: 5% initially infected, chosen 
randomly from the population. Right: 5% initially infected, chosen with probability proportional to squared degree. Each simulation curve represents 
a single simulation of the given size. As population size increases, the results converge to the theoretical prediction. 
doi:10.1371/journal.pone.0101421.g005 



solve the original equations, and then use the results to initialize 
the second set of variables. 

Case 1. We begin by reducing (3 by a factor of 2 at time t\. 
Until time t\, we are solving the original equations. By solving the 
original system until t\ we have 0rj(fi). The probability an 
individual of degree k is susceptible at time t\ is 
S{k,t l ) = S{k,0)6 0 (t l ) k . So our new !/>(*) is <J/ 1 (x)= J2k P ( k ) 
S(kfi)8 0 (t\) k x k = ^(8 0 (t\)x). We take our new d\ to have 
$l0l)=l- The intervention we are doing has no impact on the 
probability a partner is in any given state. (j>s, <f>/, and keep the 
same proportion, but are scaled up to sum to 1 so each is scaled by 
6o(h)- For example, <f> s i(?i) = 0 so ( ? l)/ e o('l) = </ , so( O ) , /''o 

(0oOi))/fo(l)0oOi). 

We restart the solutions with these new values. 

Case 2. The total probability of transmitting to a partner is 
P/ifl + j)- For this intervention we change /? so that + Y) is 
reduced by a factor of 2 at time t \ . This proceeds exactly as above 
except that the new value of ft must be smaller. 

Case 3. When we delete half the edges at random, we do not 
affect the probability that a random partner is in any given state. 
So the (j> variables rescale in the same way as for changing ft in the 
previous cases. However, \jj undergoes a more significant change. 
As a starting point, consider P\{k\), the probability an individual 
has degree k\ after edges are deleted. This depends on Po(^o), the 
probability of having kg edges prior to deletion. The relation is 



P 1 (k l )=Y l Po(ko) 



The probability the individual has degree ki and is susceptible 
at time t\ is 



«*,)= E^o)^o,O)(^)( 0 fA i ° 



So if we restart the calculations at t = t \ we have 
S(k 1 ,t l ) = Q(k 1 )/P 1 (k l ) and 



]T£Po(/co)S(A:o,0) 

A-, k n 



£Po(W/c o ,O)0o(/i) A " o £. 



k 0 \ f\\ k o- k 'i 



J2Po(ko)S(k 0 fi)6 0 (h) k o 
\l+x~ 



«Ao 0o(fi) 



(in general if we delete edges with probability p and keep with 
probability q=l—p, then the new function is \j/ \ (x) = 
\l/ 0 ([8 0 (ti)]\p + qx])). Using this new i/'ifx), the same system of 
equations holds. 

Figure 6 compares these strategies. As anticipated, the final sizes 
resulting from cases 2 and 3 are identical, regardless of the time of 
intervention (indeed this is straightforward to show from the final 
size relation). However, we see that the dynamics are significandy 
different. 

There are other ways we could capture these interventions 
mathematically. We note that in both case 1 and case 2, we could 
also capture the intervention by simply using a step change in /J. 
Case 3 could be captured by reducing (j) s , ^ and (j) R each by half 
and placing the other half into a new inactive compartment (j>x for 
the deleted edges. However, each of these is ad hoc and dependent 
on the precise details of the case. Using the approach presented 
above, we have a standard approach that will apply across a wide 
range of interventions. 

Bifurcation analysis 

We try to gain a better understanding of the epidemic threshold 
and what happens to the final size as the initial proportion infected 
is increased. Consider the final size relation found from 
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Figure 6. The impact of interventions. Epidemics begin at t=0 with 0.001 of the population infected. Left: epidemic curve without interventions, 
and with each intervention introduced at time t\ = 1.5. Right: horizontal axis is t\, showing final effectiveness if interventions introduced at different 
times 
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s(0) 



408) 

>A'(i) 



-**(0) 



+y 



P+y 



If (j) s (0) + (j) R (0) = 1, then we find that 8=1 is a solution to these 
equations. Biologically this states that if there is no infection 
initially [^/(0) = 0] there will be no infection later. To study what 
happens when (f) I (0)>0 (but possibly arbitrarily small), we begin 
by first analyzing the structure of the dynamical equation for 8 
under the assumption that <j) s (0)= 1 and <^(0) = 0, taking 6< 1. 
These assumptions contradict our initial conditions, but under- 
standing the bifurcation in this system first will lead to an easier 
understanding of the full system with <^/(0)>0. 

If (j) s (0) = \ and <j) R (0) = 0, the differential equation for 8 
becomes 



-8) 



which has 8 = 0 whenever 



MO) 
<A'(i) 



P+y 



Clearly 8=\ is an equilibrium. Close to 8=\, we write 

0=1+/*, so f(e)/ l A'(i)=i+^"(i)/>A'(i)+/' 2 iA'"(i)/2^(i)+ 

O(ji^). The value of i//"(l)/l//'(l) plays a key role in this analysis. 
Its biological interpretation comes from looking at a random 
infected individual's partner v early in the epidemic. If v is infected 
by that random infected individual, then \jl"{Y)/\j/'{\) is the 
expected number of additional partners v has (its excess degree). 
Substituting into the equation for the equilibrium we have 



\+)i=\+lill 



rm+^"'m/2+o(n 2 ) 



which yields 



fi = 0 



(5) 



So there is a bifurcation as the bracketed term passes through 
zero, when (fi-\-y)l P = \j/"(\)/\jj'(\). This is the well-known 
epidemic threshold [11]. The bifurcation is transcritical and 
corresponds to 1Zq increasing through 1. If i//"(l)/i//(l)< 
{fi + i)/P (that is TZq < 1) then fi is positive and the corresponding 
equilibrium has 6 > 1 and is unstable, while the equilibrium at 
8= 1 is stable. In our case, we will not observe 8> 1 because 6 is a 
probability. If however 4"'(X)/\jj'(X)>(fi + y)/ ft (that is TZo>l), 
then the corresponding equilibrium has 8 < 1 and is stable while 
the equilibrium at 8 = 1 is unstable. 

Previous studies found these approximate fixed points for 8, and 
used a slighdy different definition for 8 such that 0(0) was slightly 
less than 1. So the biologically implausible prediction is that for 
TZo < 1 the value of 8 increases towards the stable fixed point at 1 . 
For 1Z(, > 1 the prediction is more meaningful: 8 decreases away 
from the unstable fixed point at 1 to the lower, stable fixed point. 
When we do a more careful analysis, we now have 8(0) = 1 and 
<l>l(0)>0. The stable fixed point that was at 1 for TZq< 1 will be 
slightly decreased to 8q< 1 = 8(0). So 8 will decrease to 8$ when 
TZo < 1 • For TZq > 1 the unstable fixed point that was at 8q = 1 is 
slightly increased, so 8q > 1 = 0(0). So 8 will decrease away from 
8q. When the number of infections decays immediately (IZq < 1) 
we care about a small error in the prediction caused by the fact 
that the initial proportion infected is not exactiy zero, while if the 
number of infections grows (IZo > 1), the error caused by treating 
the initial proportion infected as asymptotically 0 is insignificant. 

We are now able to consider the effect of realistic initial 
conditions. We keep (j) R (0) = 0, but take <^/(0) to be a small positive 
number with 8(0) =1 and </> s (0)=l—(j> I (0). The bifurcation 
diagram changes slightly. Compared to the equations assuming 
<l> s (0)=l, this has the effect of decreasing 8 slightly, so the 
equilibrium values shift as shown in figure 7. In fact, stricdy 
speaking, there is no bifurcation for ^ / (0)>0. 

Below the bifurcation value of fi, the equilibrium 8 = 1 is slightly 
reduced to a 8q< 1, but remains stable. The solution with initial 
condition 8=\ converges to this equilibrium. Above the 
bifurcation, the 8=1 equilibrium is slightly increased to a 8q > 1 
and is unstable. The other equilibrium with smaller 8 is stable and 
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Figure 7. Bifurcation diagram with 0 R (O) = O, 0^(0) = 1 —0,(0), and <j>j(0) as given. The figure on right zooms in on the bifurcation point. 
Disease parameters are J?= 1 and y = 1.5. In each all members of the population have degree either 3 or 4, with the proportions chosen so that 
iJ/"(V)/iJ/(l) takes the values on the horizontal axis. The dashed curves denote unstable equilibria and the solid curves stable equilibria. Approximate 
curves (dotted) come from equation (5). Only the equilibria with /i<0 are biologically meaningful. 
doi:10.1371/journal.pone.0101421.g007 



its location is decreased slightly. Thus the solution with initial 
condition 9=1 decreases and converges to the stable solution. 

Results and Discussion 

In this paper we have extended previous mathematical methods 
for epidemic spread in networks to allow for simple models that 
can capture the impact of a non-negligible initial condition. 

Our system of equations (1) and (2) are mathematically simple 
and can be solved numerically with standard tools. Changes in the 
population's degree distribution affect but do not otherwise 
alter the structure of the equations, and in particular, the 
population can have arbitrarily large maximum degree without 
requiring any increase in the number of equations. 

Our modeling approach accurately predicts the size and 
dynamics of simulated epidemics with arbitrary sized initial 
conditions. The approach allows us to compare interventions 
introduced during the epidemic. We further see that the 
discrepancy found by [14] apparently results from the fact that 
the initial proportion infected was nonnegligible. 

We are able to investigate the dynamical structure of the 
equilibria found in the equations, showing how the epidemic 
threshold is modified when a non-negligible proportion of the 
population is infected. This helps resolve the apparent contradic- 
tion in earlier work in which 9, a measure of the probability of 
having escaped transmission, could increase in time. On biological 
grounds 9 must be monotonically decreasing. The mathematical 
inconsistency resulted from the assumption of a small initial 
condition. For IZq < 1 , in the small initial condition limit, there is 
an attracting equilibrium value where the cumulative amount of 
infection is very small. For a reduced initial infection, the 
equilibrium moves, going to no infection in the limit. The 
previous analyses effectively froze this equilibrium at its asymptotic 
limit and then considered a small, but nonzero initial amount of 
infection. So the initial condition was on the wrong side of the 
equilibrium and the mathematical model attempts to return to its 
frozen equilibrium of no cumulative infection rather than 
approaching the true equilibrium with a small number of 
recovered individuals. When TZo>^, this effect does not matter 
because the frozen equilibrium is repelling and still on the correct 
side of the initial condition. 

One of the most obvious applications of these results is to the 
understanding of the impact of an intervention that begins after a 



disease has established itself. This has been a weakness of network 
models for some time: the earliest low-dimensional models could 
only calculate static quantities such as the final size of epidemics 
assuming no intervention, while more recent approaches that 
calculate the dynamics [7,15,22] have been restricted to the 
assumption of asymptotically small initial conditions, again with no 
change in the population behavior. Because we now have a low- 
dimensional model that can account for large initial conditions, we 
can use this to restart our calculations when an intervention is to 
be implemented, or we can use the final size relation to quickly 
compare intervention effectiveness. 

We have analyzed the bifurcation structure of the final size 
relation, and used this to explain an apparent discrepancy in 
earlier work if TZq < 1 . The previous models that assumed small 
initial condition also implicitly assume that TZq > 1 . This resulted 
in a disturbing prediction for 1Zo < 1 that transmissions could be 
reversed as time progresses, and infected individuals are uninfect- 
ed. Once we correctiy account for the initial condition this 
apparent discrepancy disappears. 

If variance is large enough, then there may be a small number 
of very high degree individuals who have a macroscopic effect on 
the dynamics. Usually, increasing the population size will "drown 
out" the effects of individuals. However, when the variance is 
sufficiently large increasing the population size results in a small 
number of much higher degree individuals. These again have a 
macroscopic effect on the dynamics. Deterministic predictions will 
not be accurate: for example, how long the highest degree 
individual remains infected will influence the final size. The work 
of [8] rigorously studied the equations using small initial 
conditions, and showed that if all moments up to the fifth moment 
were finite, then these equations are accurate in the limit of a large 
network. More recently [25] has shown that the equations remain 
valid with much weaker assumptions. The equations will fail if the 
degree distribution is too broad exactly because these very rare 
individuals are able to have a population-scale impact on the 
dynamics, so the specific time and duration of their infections 
matter. 

We expect to be able to recover from this challenge however. 
After the epidemic has run for a short period of time in a 
Configuration Model network, all of these high degree individuals 
have been infected and recovered. The remaining population will 
have significantly reduced moments. At this point, the stochastic 
effects are "frozen in": the dynamics are now deterministic. We 
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can use the observed conditions at this time to initialize our new 
system of equations. 

There are some assumptions implicit in our derivation that 
deserve further attention. The model fails if 0(0), ij>s(0), <j)/(0), or 
<^(0) depend on degree of u. So if for example, we select high 
degree individuals and then infect their partners (leaving the high 
degree individuals uninfected), the model will not account for the 
fact that higher degree individuals are more likely to have infected 
partners at t = 0. The approach will fail. This is discussed in more 
detail in [26]. 

To be clear, the model does not fail if the initial individuals 
infected have higher (or lower) degree. This simply affects the 
initial conditions. Indeed, we expect that if the infection is initially 
spreading stochastically in the population, and we set t = Q to be 
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